# This code is used to estimate average delegation, constraint and discretion indices by year

library(ggplot2)
library(foreign)

data = 
  read.csv("~/Dropbox/Research/Papers/Delegation-ML-Project/FINAL-OUTPUT-DATA/APSR_RnR/discretion-delegation-constraint_MS.csv")

########### Commission data Now average by year com data

ag2<-aggregate(data, by=list(data$year), FUN=mean, na.rm=TRUE)


ones<-rep(1,dim(data)[1])

ob<-data.frame(ones,data$year)

Ns<-aggregate(ob, by = list(data$year), FUN = sum)

NLaws<-Ns$ones # This is the number of laws collected per year

## Get the standard errors of the ratios

sd.delegation<-aggregate(data$delegation,by = list(data$year), FUN = sd,na.rm=TRUE)
var.delegation<-sd.delegation$x^2
se.delegation<- sd.delegation$x/sqrt(NLaws)

sd.constraint.ms<-aggregate(data$constraint.ratio,by = list(data$year), FUN = sd,na.rm=TRUE)
var.constraint.ms<-sd.constraint.ms$x^2
se.constraint.ms<- sd.constraint.ms$x/sqrt(NLaws)

sd.discretion<-aggregate(data$discretion.index,by = list(data$year), FUN = sd,na.rm=TRUE)
var.discretion<-sd.discretion$x^2
se.discretion<- sd.discretion$x/sqrt(NLaws)


## Now let's add all of this information to the aggregated data
ag2<-data.frame(ag2,var.delegation,se.delegation,
                var.constraint.ms,se.constraint.ms,
                var.discretion,se.discretion)

EU.Laws<-NLaws
################################################################################
################################               #################################
################################  MAIN PLOTS   #################################
################################               #################################
################################################################################

setwd("~/Dropbox/Research/Papers/Delegation-ML-Project/Draft-EU-Paper/figs/APSR_RnR")

## Plot number of laws collected by  year
ggplot(ag2, aes(x=year, y=EU.Laws)) + geom_point(shape=1)+
  geom_line()+theme_classic() + ylab("EU Laws Collected") + 
  xlab("Year") +  scale_x_continuous(breaks = seq(1959, 2017, by = 2)) + 
  theme(axis.text.x = element_text(angle = 70, hjust = 1)) + geom_smooth()

ggsave("nlaws_FINAL.png")

# Restrict data to post 1963 for illustration
ag2new = ag2[ag2$year>=1970,]

## Delegation only
ggplot(ag2new, aes(x=year, y=delegation.ratio)) +
  geom_errorbar(aes(ymin=delegation.ratio-1.96*se.delegation, 
                    ymax=delegation.ratio+1.96*se.delegation),width=0.5 
                ,linetype= "dotted") +  theme_classic() + geom_point() + 
  geom_line() +
  scale_x_continuous(breaks = seq(1965, 2017, by = 2)) + 
  theme(axis.text.x = element_text(angle = 70, hjust = 1),
        legend.title = element_blank()) +
  xlab("Year") + ylab("Mean Delegation Ratio (Member States)") + geom_smooth()
#+expand_limits(y = c(0.01,0.07))

ggsave("delegation-only_FINAL_MS.png")


## Constraint only
ggplot(ag2new, aes(x=year, y=constraint.ratio)) +
  geom_errorbar(aes(ymin=constraint.ratio-1.96*se.constraint.ms, 
                    ymax=constraint.ratio+1.96*se.constraint.ms),width=0.5 
               , linetype = "dotted") +  theme_classic() + geom_point() + 
  geom_line() +
  scale_x_continuous(breaks = seq(1965, 2017, by = 2)) + 
  theme(axis.text.x = element_text(angle = 70, hjust = 1),
        legend.title = element_blank()) +
  xlab("Year") + ylab("Mean Constraint Ratio (Member States)") + geom_smooth()
#+  expand_limits(y = c(0,0.05))

ggsave("constraint-only_FINAL_MS.png")



## Discretion index
ggplot(ag2new, aes(x=year, y=discretion.index)) +
  geom_errorbar(aes(ymin=discretion.index-1.96*se.discretion, 
                    ymax=discretion.index+1.96*se.discretion),width=0.5, 
                 linetype = "dotted") +  theme_classic() + geom_point() + 
  geom_line() +
  scale_x_continuous(breaks = seq(1965, 2017, by = 2)) + 
  theme(axis.text.x = element_text(angle = 70, hjust = 1),
        legend.title = element_blank()) +
  xlab("Year") + ylab("Mean Discretion Index (Member States)")  + geom_smooth()
#+  expand_limits(y = c(0.07,0.3))

ggsave("discretion-only_FINAL_MS.png")



### reconstruct PDR and PCR
### reconstruct PDR and PCR
library(devtools)
#devtools::install_github("MatthieuStigler/RDDtools/RDDtools")
library(RDDtools)

rm(data)

data = read.csv("~/Dropbox/Research/Papers/Delegation-ML-Project/FINAL-OUTPUT-DATA/APSR_RnR/final_delegation_data_rnr.csv")

# Need to create a pdr plot for each constraint category
#[1] "X.1"                                "X"                                  "newprovision.celex"                
#[4] "newprovisiondata.month"             "newprovision.year"                  "ms.appeals.procedures.probs"       
#[7] "ms.appeals.procedures.preds"        "ms.consultation.requirements.probs" "ms.consultation.requirements.preds"
#[10] "ms.executive.action.possible.probs" "ms.executive.action.possible.preds" "ms.executive.action.required.probs"
#[13] "ms.executive.action.required.preds" "ms.reporting.requirements.probs"    "ms.reporting.requirements.preds"   
#[16] "ms.rulemaking.requirements.probs"   "ms.rulemaking.requirements.preds"   "ms.spending.limits.probs"          
#[19] "ms.spending.limits.preds"           "ms.time.limit.probs"                "ms.time.limit.preds"               
#[22] "delegation.probs"                   "delegation.predictions"   
attach(data)



appeals.index = which(ms.appeals.procedures.preds == 1)
consultation.index = which(ms.consultation.requirements.preds == 1)
executive.action.possible.index = which(ms.executive.action.possible.preds == 1)
executive.action.req.index = which(ms.executive.action.required.preds == 1)
reporting.index = which(ms.reporting.requirements.preds == 1)
rulemaking.index = which(ms.rulemaking.requirements.preds == 1)
spending.index = which(ms.spending.limits.preds == 1)
timelimit.index = which(ms.time.limit.preds == 1)
delegation.index = which(delegation.predictions == 1)


appeals.pcr = ms.appeals.procedures.probs[appeals.index]
consultation.pcr = ms.consultation.requirements.probs[consultation.index]
executive.action.pcr = ms.executive.action.possible.probs[executive.action.index]
executive.action.req.pcr = ms.executive.action.required.probs[executive.action.req.index]
reporting.pcr = ms.reporting.requirements.probs[reporting.index]
rulemaking.pcr = ms.rulemaking.requirements.probs[rulemaking.index]
spending.pcr = ms.spending.limits.probs[spending.index]
timelimit.pcr = ms.time.limit.probs[timelimit.index]
delegation.pdr = delegation.probs[delegation.index]

#### Now let's estimate a bunch of models
y1 = newprovision.year[appeals.index]
y2 = newprovision.year[consultation.index]
y3 = newprovision.year[executive.action.possible.index]
y4 = newprovision.year[executive.action.req.index]
y5 = newprovision.year[reporting.index]
y6 = newprovision.year[rulemaking.index]
y7 = newprovision.year[spending.index]
y8 = newprovision.year[timelimit.index]
y9 = newprovision.year[delegation.index]


d1 = RDDdata(y = appeals.pcr, x = y1, cutpoint = 1993.5)
d2 = RDDdata(y = consultation.pcr, x =y2, cutpoint = 1993.5)
d3 = RDDdata(y = executive.action.pcr, x =y3, cutpoint = 1993.5)
d4 = RDDdata(y = executive.action.req.pcr, x =y4, cutpoint = 1993.5)
d5 = RDDdata(y = reporting.pcr, x =y5, cutpoint = 1993.5)
d6 = RDDdata(y = rulemaking.pcr, x =y6, cutpoint = 1993.5)
d7 = RDDdata(y = spending.pcr, x =y7, cutpoint = 1993.5)
d8 = RDDdata(y = timelimit.pcr, x =y8, cutpoint = 1993.5)
d9 = RDDdata(y = delegation.pdr, x = y9 , cutpoint = 1993.5)

bw_ik1 = RDDbw_IK(d1)
bw_ik2 = RDDbw_IK(d2)
bw_ik3 = RDDbw_IK(d3)
bw_ik4 = RDDbw_IK(d4)
bw_ik5 = RDDbw_IK(d5)
bw_ik6 = RDDbw_IK(d6)
bw_ik7 = RDDbw_IK(d7)
bw_ik8 = RDDbw_IK(d8)
bw_ik9 = RDDbw_IK(d9)

r1 = RDDreg_lm(RDDobject = d1, bw = bw_ik1,order=3)
r2 = RDDreg_lm(RDDobject = d2, bw = bw_ik2,order=3)
r3 = RDDreg_lm(RDDobject = d3, bw = bw_ik3,order=3)
r4 = RDDreg_lm(RDDobject = d4, bw = bw_ik4,order=3)
r5 = RDDreg_lm(RDDobject = d5, bw = bw_ik5,order=3)
r6 = RDDreg_lm(RDDobject = d6, bw = bw_ik6,order=3)
r7 = RDDreg_lm(RDDobject = d7, bw = bw_ik7,order=3)
r8 = RDDreg_lm(RDDobject = d8, bw = bw_ik8,order=3)
r9 = RDDreg_lm(RDDobject = d9, bw = bw_ik9,order=3)

## Finally let's plot all of them and save them

setwd("~/Dropbox/Research/Papers/Delegation-ML-Project/Draft-EU-Paper/APSR_RnR/figs")

png("pcr1_FINAL_ms.png")
plot(r1, xlab = "Year", ylab = "Mean Constraint Probability (Appeals)")
abline(v = 1993.5,lty="dotted")
dev.off()

png("pcr2_FINAL_ms.png")
plot(r2, xlab = "Year", ylab = "Mean Constraint Probability (Consultation)")
abline(v = 1993.5,lty="dotted")
dev.off()

png("pcr3_FINAL_ms.png")
plot(r3, xlab = "Year", ylab = "Probabilistic Constraint Ratio (Executive Action Possible))")
abline(v = 1993.5,lty="dotted")
dev.off()

png("pcr4_FINAL_ms.png")
plot(r4, xlab = "Year", ylab = "Mean Constraint Probability (Executive Action Required)")
abline(v = 1993.5,lty="dotted")
dev.off()

png("pcr5_FINAL_ms.png")
plot(r5, xlab = "Year", ylab = "Mean Constraint Probability (Reporting)")
abline(v = 1993.5,lty="dotted")
dev.off()

png("pcr6_FINAL_ms.png")
plot(r6, xlab = "Year", ylab = "Mean Constraint Probability (Rulemaking)")
abline(v = 1993.5,lty="dotted")
dev.off()

png("pcr7_FINAL_ms.png")
plot(r7, xlab = "Year", ylab = "Mean Constraint Probability (Spending)")
abline(v = 1993.5,lty="dotted")
dev.off()

png("pcr8_FINAL_ms.png")
plot(r8, xlab = "Year", ylab = "Mean Constraint Probability (Time Limits)")
abline(v = 1993.5,lty="dotted")
dev.off()

png("pdr_FINAL_ms.png")
plot(r9, xlab = "Year", ylab = "Mean Delegation Probability (Authority)",ylim = c(0.5,1))
abline(v = 1993.5,lty="dotted")
dev.off()


####################################################################################################
####################################################################################################
####################################################################################################
################################Amending v Non-Amending#############################################
####################################################################################################
####################################################################################################
####################################################################################################

#appeals.pcr = ms.appeals.procedures.probs[appeals.index]
#consultation.pcr = ms.consultation.requirements.probs[consultation.index]
#executive.action.pcr = ms.executive.action.possible.probs[executive.action.index]
#executive.action.req.pcr = ms.executive.action.required.probs[executive.action.req.index]
#reporting.pcr = ms.reporting.requirements.probs[reporting.index]
#rulemaking.pcr = ms.rulemaking.requirements.probs[rulemaking.index]
#spending.pcr = ms.spending.limits.probs[spending.index]
#timelimit.pcr = ms.time.limit.probs[timelimit.index]
#delegation.pdr = delegation.probs[delegation.index]

data = read.csv("~/Dropbox/Research/Papers/Delegation-ML-Project/FINAL-OUTPUT-DATA/APSR_RnR/final_delegation_data_rnr.csv")

# Need to create a pdr plot for each constraint category
#[1] "X.1"                                "X"                                  "newprovision.celex"                
#[4] "newprovisiondata.month"             "newprovision.year"                  "ms.appeals.procedures.probs"       
#[7] "ms.appeals.procedures.preds"        "ms.consultation.requirements.probs" "ms.consultation.requirements.preds"
#[10] "ms.executive.action.possible.probs" "ms.executive.action.possible.preds" "ms.executive.action.required.probs"
#[13] "ms.executive.action.required.preds" "ms.reporting.requirements.probs"    "ms.reporting.requirements.preds"   
#[16] "ms.rulemaking.requirements.probs"   "ms.rulemaking.requirements.preds"   "ms.spending.limits.probs"          
#[19] "ms.spending.limits.preds"           "ms.time.limit.probs"                "ms.time.limit.preds"               
#[22] "delegation.probs"                   "delegation.predictions"   
attach(data)



appeals.index = which(ms.appeals.procedures.preds == 1)
consultation.index = which(ms.consultation.requirements.preds == 1)
executive.action.possible.index = which(ms.executive.action.possible.preds == 1)
executive.action.req.index = which(ms.executive.action.required.preds == 1)
reporting.index = which(ms.reporting.requirements.preds == 1)
rulemaking.index = which(ms.rulemaking.requirements.preds == 1)
spending.index = which(ms.spending.limits.preds == 1)
timelimit.index = which(ms.time.limit.preds == 1)
delegation.index = which(delegation.predictions == 1)


appeals.pcr = ms.appeals.procedures.probs[appeals.index]
consultation.pcr = ms.consultation.requirements.probs[consultation.index]
executive.action.pcr = ms.executive.action.possible.probs[executive.action.index]
executive.action.req.pcr = ms.executive.action.required.probs[executive.action.req.index]
reporting.pcr = ms.reporting.requirements.probs[reporting.index]
rulemaking.pcr = ms.rulemaking.requirements.probs[rulemaking.index]
spending.pcr = ms.spending.limits.probs[spending.index]
timelimit.pcr = ms.time.limit.probs[timelimit.index]
delegation.pdr = delegation.probs[delegation.index]



# Compare amending probabilities for each categority with non-amending.
amendingdata = read.csv("~/Dropbox/Research/Papers/Delegation-ML-Project/FINAL-OUTPUT-DATA/AmendingID/amending-id-celex.csv")

amendingid= amendingdata$amending.celex

rm(amendingdata)

amendingdummy = rep(0,length(newprovision.celex))

nrows = length(newprovision.celex)

for(i in 1:nrows){
  print(i)
  if(sum(newprovision.celex[i] ==amendingid ) ==1){
    amendingdummy[i] <- 1
    amendingid = amendingid[-c(i)] # Remove the id for faster computation
  }
}

# Write out a new dataframe for the future: 
final_delegation_data_rnr_amending = data.frame(data,amendingdummy)
write.csv(final_delegation_data_rnr_amending, "~/Dropbox/Research/Papers/Delegation-ML-Project/FINAL-OUTPUT-DATA/APSR_RnR/final_delegation_data_rnr_amending.csv" )
rm(final_delegation_data_rnr_amending)

appeals.index.amending = which(ms.appeals.procedures.preds == 1&amendingdummy == 1)
consultation.index.amending = which(ms.consultation.requirements.preds == 1&amendingdummy == 1)
executive.action.possible.index.amending = which(ms.executive.action.possible.preds == 1&amendingdummy == 1)
executive.action.req.index.amending = which(ms.executive.action.required.preds == 1&amendingdummy == 1)
reporting.index.amending = which(ms.reporting.requirements.preds == 1&amendingdummy == 1)
rulemaking.index.amending = which(ms.rulemaking.requirements.preds == 1&amendingdummy == 1)
spending.index.amending = which(ms.spending.limits.preds == 1&amendingdummy == 1)
timelimit.index.amending = which(ms.time.limit.preds == 1&amendingdummy == 1)
delegation.index.amending = which(delegation.predictions == 1&amendingdummy == 1)


appeals.index.noamending = which(ms.appeals.procedures.preds == 1&amendingdummy == 0)
consultation.index.noamending = which(ms.consultation.requirements.preds == 1&amendingdummy == 0)
executive.action.possible.index.noamending = which(ms.executive.action.possible.preds == 1&amendingdummy == 0)
executive.action.req.index.noamending = which(ms.executive.action.required.preds == 1&amendingdummy == 0)
reporting.index.noamending = which(ms.reporting.requirements.preds == 1&amendingdummy == 0)
rulemaking.index.noamending = which(ms.rulemaking.requirements.preds == 1&amendingdummy == 0)
spending.index.noamending = which(ms.spending.limits.preds == 1&amendingdummy == 0)
timelimit.index.noamending = which(ms.time.limit.preds == 1&amendingdummy == 0)
delegation.index.noamending = which(delegation.predictions == 1&amendingdummy == 0)

# Create amending v non-amending probabilities
#amending
appeals.pcr.amending = ms.appeals.procedures.probs[appeals.index.amending]
consultation.pcr.amending = ms.consultation.requirements.probs[consultation.index.amending]
executive.action.pcr.amending = ms.executive.action.possible.probs[executive.action.index.amending]
executive.action.req.pcr.amending = ms.executive.action.required.probs[executive.action.req.index.amending]
reporting.pcr.amending = ms.reporting.requirements.probs[reporting.index.amending]
rulemaking.pcr.amending = ms.rulemaking.requirements.probs[rulemaking.index.amending]
spending.pcr.amending = ms.spending.limits.probs[spending.index.amending]
timelimit.pcr.amending = ms.time.limit.probs[timelimit.index.amending]
delegation.pdr.amending = delegation.probs[delegation.index.amending]

#no amending
appeals.pcr.noamending = ms.appeals.procedures.probs[appeals.index.noamending]
consultation.pcr.noamending = ms.consultation.requirements.probs[consultation.index.noamending]
executive.action.pcr.noamending = ms.executive.action.possible.probs[executive.action.index.noamending]
executive.action.req.pcr.noamending = ms.executive.action.required.probs[executive.action.req.index.noamending]
reporting.pcr.noamending = ms.reporting.requirements.probs[reporting.index.noamending]
rulemaking.pcr.noamending = ms.rulemaking.requirements.probs[rulemaking.index.noamending]
spending.pcr.noamending = ms.spending.limits.probs[spending.index.noamending]
timelimit.pcr.noamending = ms.time.limit.probs[timelimit.index.noamending]
delegation.pdr.noamending = delegation.probs[delegation.index.noamending]

t1 = t.test(delegation.pdr.amending,delegation.pdr.noamending)
t2 = t.test(timelimit.pcr.amending,timelimit.pcr.noamending)
t3 = t.test(spending.pcr.amending,spending.pcr.noamending)
t4 = t.test(rulemaking.pcr.amending,rulemaking.pcr.noamending)
t5 = t.test(reporting.pcr.amending,reporting.pcr.noamending)
t6 = t.test(executive.action.req.pcr.amending,executive.action.req.pcr.noamending)
t8 = t.test(consultation.pcr.amending,consultation.pcr.noamending)
t9 = t.test(appeals.pcr.amending,appeals.pcr.noamending)



dfreedom = c(t1$parameter,
            t2$parameter,
            t3$parameter,
            t4$parameter,
            t5$parameter,
            t6$parameter,
            t8$parameter,
            t9$parameter)

dfreedom = as.vector(round(dfreedom,digits = 0))

amendingest = c(t1$estimate[1],
                t2$estimate[1],
                t3$estimate[1],
                t4$estimate[1],
                t5$estimate[1],
                t6$estimate[1],
                t8$estimate[1],
                t9$estimate[1])

noamendingest = c(t1$estimate[2],
                t2$estimate[2],
                t3$estimate[2],
                t4$estimate[2],
                t5$estimate[2],
                t6$estimate[2],
                t8$estimate[2],
                t9$estimate[2])

diff = noamendingest-amendingest

pvalues = c(t1$p.value,
                  t2$p.value,
                  t3$p.value,
                  t4$p.value,
                  t5$p.value,
                  t6$p.value,
                  t8$p.value,
                  t9$p.value)






library(xtable)

amendingcompare = data.frame(noamendingest,amendingest,diff,pvalues,dfreedom)

names(amendingcompare) = c("Non-Amending", "Amending", "Difference", "P Value","Degrees of Freedom")

rownames(amendingcompare) = c("Mean Pr(Delegation)",
             "Mean Pr(Time Limit Constraint)",
             "Mean Pr(Spending Constraint)",
             "Mean Pr(Rulemaking Constraint)",
             "Mean Pr(Reporting Constraint)",
             "Mean Pr(Exec. Action Req Constraint)",
             "Mean Pr(Cons. Constraint)",
             "Mean Pr(Appeals Constraint)"
)

digits(amendingcompare)<-c(0,3,3,3,3,0)

xtable(amendingcompare,digits = c(0,3,3,3,3,0))


## Plots of predicted constraint v actual in Franchino
franchino.ratio<-
  read.csv("~/Dropbox/Research/Papers/Delegation-ML-Project/Data/delegationfran.csv")


























